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Jamming is a geometric phase transition occurring in dense particle systems in the absence of 
temperature. We use computer simulations to analyse the effect of thermal fluctuations on several 
signatures of the transition. We show that scaling laws for bulk and shear moduli only become rele¬ 
vant when thermal fluctuations are extremely small, and propose their relative ratio as a quantitative 
signature of jamming criticality. Despite the nonequilibrium nature of the transition, we find that 
thermally induced fluctuations and mechanical responses obey equilibrium fluctuation-dissipation re¬ 
lations near jamming, provided the appropriate fluctuating component of the particle displacements 
is analysed. This shows that mechanical moduli can be directly measured from particle positions 
in mechanically unperturbed packings, and suggests that the definition of a “nonequilibrium index” 
is unnecessary for amorphous materials. We find that fluctuations of particle displacements are 
spatially correlated, and define a transverse and a longitudinal correlation lengthscales which both 
diverge as the jamming transition is approached. We analyse the frozen component of density fluc¬ 
tuations and find that it displays signatures of nearly-hyperuniform behaviour at large lengthscales. 

This demonstrates that hyperuniformity in jammed packings is unrelated to a vanishing compress¬ 
ibility and explains why it appears remarkably robust against temperature and density variations. 

Differently from jamming criticality, obstacles preventing the observation of hyperuniformity in 
colloidal systems do not originate from thermal fluctuations. 

PACS numbers: 05.10.-a, 05.20.Jj, 64.70.qj 


I. INTRODUCTION 

A jamming transition M occurs when it becomes 
too difficult to compress any further a dense assembly of 
hard objects, whose compressibility then vanishes. Re¬ 
markably, the same transition also controls the loss of 
mechanical rigidity observed when soft particles are de¬ 
compressed in the absence of thermal fluctuations, such 
as foams or emulsion droplets, whose bulk modulus then 
vanishes. In both cases, the key variable controlling the 
physics is the particle connectivity, the jamming transi¬ 
tion corresponding to the isostatic situation where just 
enough contacts are present to insure mechanical stabil¬ 
ity The distance to isostaticity controls the diver¬ 

gence of mechanical moduli observed when compressing 
hard particles and their vanishing when decompressing 
soft particles 3 ■ The deep connection between geom¬ 
etry and mechanical responses shows that the criticality 
observed in the vicinity of the transition follows from the 
unambiguous identification of particle contacts. There¬ 
fore, when thermal fluctuations are present, for instance 
when considering colloidal particles (such as microgels, 
emulsions, PMMA colloids), the contacts can be ‘blurred’ 
by thermal agitation and cannot be resolved, which chal¬ 
lenges the possibility to observe jamming criticality in 
experiments. In other words, a thermal system cannot 
know how far it is from isostaticity , and the associated 
criticality is easily destroyed by temperature [la¬ 
in previous work, the role of thermal fluctuations near 
jamming has been explored to understand the influence of 
finite temperatures on various physical quantities such as 
microscopic dynamics, microstructure, contact number, 


mechanical properties [IE- In particular, a computer 
study of the single particle dynamics revealed the exis¬ 
tence of a very narrow region in the (density, tempera¬ 
ture) phase diagram where jamming criticality can be ob¬ 
served, which excludes most colloidal studies to date El 
. More recent experiments have concentrated on collec¬ 
tive static properties, such as mechanical shear and bulk 
moduli and structure factors, and the results were anal¬ 
ysed using power laws that are valid, strictly speaking, for 
fully athermal systems (U] . To assess the validity of this 
description, one needs to extend the analysis of Ref. El 
to mechanical moduli to understand whether their criti¬ 
cal behaviour is robust against thermal fluctuations. The 
first goal of our work is to analyse the effect of thermal 
fluctuations on mechanical moduli near jamming. 

Another property characterizing jammed packings is 
their hyperuniformity, which was revealed by analysing 
the large-distance scaling of volume fraction fluctuations 
E}. For monodisperse spherical particles of diameter 
cr, this reduces to studying the ordinary static structure 
factor, S(k), whose low-wavevector behavior obeys a non¬ 
trivial, characteristic linear behaviour, S(ka -C 1) ~ k, 
which shows that density fluctuations are suppressed at 
large scale E, El- This behavior has been observed 
numerically in particle packings prepared exactly at the 
jamming transition [22H261I , and experimentally in ather¬ 
mal granular materials [25[. Experiments performed with 
colloidal particles appear challenging and report only 
very weak signs of hyperuniformity j27H29|- A possi¬ 
ble explanation could be that hyperuniform behaviour is 
blurred by thermal fluctuations, as are other signatures 
of the jamming transition. However, hyperuniformity is 
a property of the packings at large lengthscale and the 


2 


above argument regarding the resolution of particle con¬ 
tacts is not obviously relevant. Therefore, if hyperunifor¬ 
mity were affected by thermal fluctuations acting at the 
(vanishingly small) contact lengthscale, it would directly 
establish that hyperuniformity is another critical prop¬ 
erty associated to the jamming transition. The second 
goal of our work is to test whether hyperuniformity is ro¬ 
bust against thermal fluctuations, and, more fundamen¬ 
tally, whether hyperuniformity is deeply related to the 
jamming transition, or is instead a distinct phenomenon. 

Thermal fluctuations in jammed packings not only 
raise practical issues about experimental observations, 
they also pose fundamental challenges related to the 
nonequilibrium nature of the jamming transition. As 
mentioned above, mechanical moduli display power law 
behaviour near jamming at zero temperature. How¬ 
ever, for materials at thermal equilibrium, mechanical 
response functions are directly related to mechanical fluc¬ 
tuations induced by thermal motion and thus to equilib¬ 
rium structure factors through fluctuation-dissipation re¬ 
lations BIU]. Near the nonequilibrium jamming tran¬ 
sition at finite temperatures, two behaviors are then pos¬ 
sible: 

1) Fluctuations and responses do not obey equilibrium 
relations , so that mechanical moduli and structure fac¬ 
tors have independent density and temperature depen¬ 
dences. This hypothesis suggests that it could be use¬ 
ful to introduce novel variables to quantify deviations 
from e quili brium relations, such as nonequilibrium in¬ 
dex [32|, |33( or effective temperatures [3J, |35| , generically 
defined as ratios between fluctuations and responses. 
In that case, structure factors live an independent life 
from mechanical responses, and they may display an 
independent set of critical properties, but they may 
also have unremarkable behaviour near jamming. This 
general hypothesis has been advocated in particular in 
Refs. [ 32 I [331, where a diverging nonequilibrium index 
and a diverging nonequilibrium lengthscale were defined 
from fluctuations and responses of hard sphere systems 
approaching jam ming [32j], and later extended to generic 
amorphous solids I33j . 

2) Fluctuations and responses obey equilibrium rela¬ 
tions, and the critical behaviour of mechanical moduli 
should have a counterpart in fluctuations of particle po¬ 
sitions and structure factors. In that case, if thermal fluc¬ 
tuations are finite (but still sufficiently small that they 
do not blur the jamming criticality!), interesting critical 
behaviour should be observed in density fluctuations of 
mechanically unperturbed packings. In particular, one 
may expect the emergence of diverging lengthscales in 
collective structure factors of jammed materials. 

The third goal of our work is to decide which of the two 
above scenarios is valid, and whether interesting length- 
scales and nonequilibrium indicators emerge from the 
analysis of structure factors. 

To achieve our three main goals, we use computer sim¬ 
ulations of a simple model of soft harmonic particles §56} 
to analyse the influence of thermal fluctuations on me¬ 


chanical responses and structure factors in the vicinity 
of the jamming transition. Harmonic spheres are conve¬ 
nient because they allow studies of multiple, experimen¬ 
tally relevant, routes to jam ming in the (density, temper¬ 
ature) phase diagram [ 3 , [lj, |36 h 39} . Therefore, a single 
model provides us with decisive answers to the three sets 
of questions mentioned above, that can be summarized 
as follows. 

(i) We find that mechanical moduli are as sensitive to 
thermal fluctuations as single particle dynamics and their 
associated power law behaviour is not a good starting 
point to theoretically describe existing colloidal experi¬ 
ments. 

(ii) By contrast, hyperuniformity is extremely robust 
to the addition of thermal perturbations, and even to 
changes in packing fraction, suggesting that it should 
in fact be far easier to observe in experiments than the 
jamming criticality, even though the present state of the 
literature suggests the opposite. We also conclude that 
hyperuniformity bears no deep relation to the jamming 
transition, and in particular we show that it is fully un¬ 
related to the critical behavior of the mechanical com¬ 
pressibility. 

(iii) Equilibrium fluctuation-dissipation relations are 
perfectly obeyed near jamming, suggesting it is unnec¬ 
essary to define quantitative indicators for the degree of 
‘nonequilibriunmess’ near jamming. It also implies that 
structure factors display critical properties and reveal di¬ 
verging lengthscales, that we define, analyse, and com¬ 
pare to previously studied critical lengthscales. 

This article is organised as follows. In Sec. [IT] we define 
the model and our numerical strategy. In Sec. nm we 
analyse the behaviour of mechanical moduli. In Sec. IIVI 
we define and analyse the behavior of structure factors 
and their associated lengthscales. In Sec. |V] we discuss 
the hyperuniformity of jammed packings. In Sec. IVII we 
summarize and discuss our results. 


II. MODEL AND SIMULATION 

We consider a system of monodisperse harmonic 
spheres, interacting through a pairwise potential §56} , 

v(nj) = |(1 - rij/a) 2 e(a - nj), ( 1 ) 

where Q(x) is the Heaviside function, r l3 is the distance 
between particles i and j, and a is the particle diameter. 
Throughout this work, length, energy, temperature and 
mechanical moduli are measured in units of a, e, e/ks , 
and e/a 3 , respectively. 

We use molecular dynamics simulations [40j | to com¬ 
pute the mechanical moduli and static structure fac¬ 
tors of the system at finite temperature. The setting 
of the calculations is essentially similar to our previous 
work [11]. We first generate a random configuration of 
N = 64,000 particles in a simulation box with peri¬ 
odic boundary conditions. The linear dimension of the 


3 


box L is adjusted to realize the packing fraction frac¬ 
tion p = 0.80, where p = na 3 N/(6L 3 ). Starting from 
this configuration, we perforin molecular dynamics simu¬ 
lations at T = 10 — 5 , where we integrate Newton’s equa¬ 
tions of motion using velocity rescaling to control the 
temperature. This can be seen as an extensive aging of 
the system starting from T = oo down to T = 10 -5 . We 
find that temperature is low enough that aging dynam¬ 
ics eventually stops and the energy and average particle 
positions reach well-defined values that do not depend 
on waiting time any longer. After this long annealing of 
the system, we change the density and temperature to 
the desired values smoothly, letting the system relax at 
each state point before taking any measurement. This 
protocol allows us to study essentially the same particle 
packing at different densities and temperatures. We do 
not study temperatures larger than T = 10” 5 , and in the 
studied regime particle diffusion and rearrangements can 
be safely neglected. 

At each state point, we then perform molecular dy¬ 
namics simulations to calculate the mechanical moduli 
and static structure factors, where we integrate Newton’s 
equations of motion in the NVE ensemble, i.e. without 
thermostat. We denote the long-time average in these 
calculations with brackets, (••■). For the particular con¬ 
figuration which is analysed extensively in this work, the 
jamming density is pj — 0.648. 

For the specific purposes of Sec. El we also calculate 
the static structure factor of jammed harmonic spheres 
with a larger number of particles, N = 512,000, at 
strictly zero temperature. To this end, we generate a 
random configuration of particles in a simulation box 
with p = 0.80, and then apply the FIRE algorithm to 
minimize the potential energy of the system at this den¬ 
sity HU. Starting from this jammed configuration, we 
decrease the density by small steps and minimize the 
potential energy after each step to obtain a series of 
jammed particle configurations over a range of densities 
Q. To increase the statistics of the results obtained for 
these zero-temperature packings, we followed this pro¬ 
cedure starting from 8 independent random configura¬ 
tions, and finally averaged the results over these inde¬ 
pendent runs. For this series of simulations, we find that 
pj ~ 0.64571 ±0.00012, where the errorbar indicates the 
standard deviation among independent packings. Av¬ 
eraging over those configurations is therefore accurate as 
long as the distance to the jamming density is larger than 
\ip — tpj\ ss 0.0001. 


III. MECHANICAL MODULI AND JAMMING 
CRITICALITY 


A. Bulk modulus 

We start our analysis with the calculation of the bulk 
modulus. The isothermal bulk modulus B quantifies the 
resistance of the system to compression. Its definition 
then naturally involves the pressure (noted P) derivative 
of the volume (noted V), 

B - V (w)r < 2 > 

We first calculate B through the response formula 
Eq. © , where the pressure P is calculated from the virial 
formula 


P = pT + 


{W) 
V ’ 


(3) 


where W = Yij r ij v '{ r ij)/3 is the virial |31J. The bulk 
modulus is of course inversely proportional to the isother¬ 
mal compressibility, B = 1 /\ t - In practice we mea¬ 
sure the pressure for various densities, and estimate the 
derivative in Eq. © using finite differences, which sug¬ 
gests that compressions or decompressions yield the same 
results. The numerical results are shown with open sym¬ 
bols and dashed lines in Fig. [T| 

The density dependence of the bulk modulus strongly 
depends on the temperature. At lower temperature, e.g. 
T = 10 , the bulk modulus increases very sharply with 

density when ip < pj, and becomes essentially density- 
independent when p > ipj. This behavior can be un¬ 
derstood as a smooth crossover between the jamming 
of Brownian hard spheres and the unjamming of non- 
Brownian soft spheres. For p < ipj and T —> 0, the 
particles have vanishing overlaps and essentially explore 
hard sphere configurations. The pressure of hard spheres 
diverges at the jamming transition as P ~ T(pj — p) _1 , 
and as a result the bulk modulus behaves as (42[ 


B ~ T(pj - p)~ 2 , (4) 


which underlies both the critical nature of the transition 
and the entropic origin of solidity in hard particle sys¬ 
tems. Our numerical results for p < pj can be well fit¬ 
ted with this power-law divergence. Equivalently, Eq. ® 
implies that the isothermal compressibility Xt vanishes 
quadratically with (pj — p) in this regime. 

On the other hand, for p > pj at lower temperatures, 
particles may have finite overlaps, and thermal fluctua¬ 
tions play a small role. Thus, the system corresponds to 
non-Brownian soft spheres. For this system, the pres¬ 
sure emerges continuously at the jamming transition, 
P ~ Up — pj), and thus the bulk modulus is expected 
to be (7] 


In this section we analyse the temperature and den¬ 
sity dependences of bulk and shear moduli of harmonic 
spheres in the vicinity of the jamming transition occur¬ 
ring at (T = 0 , p = pj). 


B ~ const. (5) 

Again, this behavior is in good agreement with the results 
shown in Fig. [Q 
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FIG. 1: (Color online) Density dependence of various quanti¬ 
ties for different temperatures indicated in the label, (a) Bulk 
modulus. Results from the pressure derivative, Eq. m , are 
shown with open symbols and dashed lines, results from the 
fluctuation formula, Eq. are shown with filled symbols 
and solid lines, (b) Shear modulus obtained from Eq. 0. 
Open diamonds indicate data for T = 0 obtained from re¬ 
sponse functions [?:]. (c) Ratio of the bulk to shear modulus. 
The vertical dashed line indicates the location of the T = 0 
jamming transition at ipj ~ 0.648. 


When temperature is increased above T = 10 -8 , the 
difference of behaviour observed on both sides of the tran¬ 
sitions becomes smaller, and for the largest studied tem¬ 
perature, T = 10 -5 , the bulk modulus is a smooth func¬ 
tion of density across ipj. This shows that systems char¬ 
acterized by T = 10 -5 (in units of the particle softness e 
0 ) are unable to reveal signs of the underlying jamming 


transition at T = 0. In particular, the solid behaviour 
of these systems is better interpreted as resulting from 
hitting a glass transition line T g (ip) [43j. This qualitative 
behavior is consistent with our previous discussion of the 
single particle dynamics & 

We now use a different approach to compute the bulk 
modulus which does not involve a response function, but 
stems instead from studying thermal fluctuations in a me¬ 
chanically unperturbed material. Within the framework 
of equilibrium statistical mechanics, the bulk modulus 
can be directly related to the fluctuations of the pres¬ 
sure. In the NVE ensemble that we use to compute the 
pressure fluctuations, the formula for the isothermal bulk 
modulus reads [40j 


B = P 


<W2> <P 2 } - <P) 2 


V 


T 


■V- 


V-M, 

6 pCy 


( 6 ) 


where W 2 = + r‘- j v"(r ij ))/9 is the hyper- 

virial, cy is the specific heat per particle, and 7 y is the 
thermal pressure coefficient [40j. The last term in Eq. ([ 6 ]) 
arises because we work in the microcanonical ensemble 
where the energy is conserved. 

We have calculated the bulk modulus through Eq. ©, 
and report the results as filled symbols in Fig. [I] Clearly, 
the results are in excellent agreement with the ones ob¬ 
tained from the response formula, Eq. ©■ Therefore, we 
conclude that fluctuations and response functions yield 
identical results for repulsive colloidal particles near jam¬ 
ming over a broad range of densities and temperatures. 
The natural interpretation is that equilibrium relations 
appear satisfied because deeply jammed solids dynami¬ 
cally explore a restricted portion of the configurational 
space located near a metastable amorphous configura¬ 
tion. In other words, the system is locally in equilib¬ 
rium, even though ergodicity is globally broken. Using 
the lan gua ge of the two-temperature scenario for aging 
glasses [3J|, the thermal fluctuations that are probed in 
the present system correspond to the fast degrees of free¬ 
dom that appear locally equilibrated at the temperature 
of the thermal bath. This physical perspective justifies 
why it is unnecessary to introduce an effective tempera¬ 
ture for the slow degrees of freedom, because these are 
completely frozen in the type of analysis that we per¬ 
form. In other words, vibrational motion does not reveal 
the nonequilibrium nature of the glass. 

Our results seem to contradict previous work f32| in¬ 
troducing a nonequilibrium index X to quantify devia¬ 
tion from equilibrium behaviour between the bulk mod¬ 
ulus B (measured as a reponse function) and the small- 
wavevector limit of the static structure factor, S(k —> 0), 
which reduces to the isothermal compressibility pT\T 
for equilibrium fluids. In Sec. EY1 we clarify the re¬ 
lation between structure factors, pressure fluctuations, 
and compressibility, and show that the physical content 


of the “nonequilibrium” index introduced in Ref. 32] 


can in fact be fully understood in terms of equilibrium 
fluctuation-dissipation relations. 
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B. Shear modulus 

We now turn to the analysis of the shear modulus, 
G. There are several ways to compute the shear modu¬ 
lus in numerical simulations (40]. The first option is to 
use global fluctuations. Just as the bulk modulus can 
be determined from the fluctuations of the pressure, the 
shear modulus can be obtained from the fluctuations of 
the shear stress. We have first tried to use this approach 
to calculate the shear modulus, but found that an ac¬ 
curate determination of G is not easy because the fluc¬ 
tuation formula requires to take the difference between 
large numbers which largely cancel and have important 
statistical fluctuations. 

To overcome this problem, we use the alternative 
method introduced in Ref. (44|. In this approach, the 
shear modulus is calculated as the k —> 0 limit of the 
the correlation function of the transverse displacement 
S T (k), 


G = 


r P T 
Lo St{ k)‘ 


(7) 


The precise definition and detailed analysis of Sr{k) will 
be given in Sec. IIV1 see Eq. m- For the moment, we 
simply notice that the determination of G from Eq. 0 
clearly stems from spontaneous fluctuations, and this ap¬ 
proach thus differs from earlier determinations based on 
response functions |2]. Here, we concentrate on the tem¬ 
perature and density dependences of the shear modulus 
G, and report our results in Fig. [L] Note that this defi¬ 
nition of the shear modulus does not require testing the 
validity of linear response, and does not depend either 
on the chosen direction for shearing. Although the shear 
modulus measured as a response function may depend on 
the direction of shear j4f| , all the directions of the shear 
modes are averaged out in the definition of Sr(k) that 
we use in Eq. (Ell- 

First, we check the validity of the fluctuation formula 
Eq. 0 . To this end we compare our results to the shear 
modulus obtained from the response function in Fig. |T| 
Although the available data is limited to the density 
above ipj at T = 0 Q, our results at lower tempera¬ 
ture are quantitatively the same as the data from the 
response function at T = 0. This confirms that fluctu¬ 
ations and response functions yield identical results for 
the shear modulus as well. A similar agreement between 
response and correlations for the shear modulus was re¬ 
ported in other glassy systems [HI, HH , which appears as 
a robust result. 

The overall behavior of the shear modulus is quali¬ 
tatively similar to the one of the bulk modulus. At 
lower temperature, the shear modulus also increases very 
sharply with <p below the jamming density, and has a 
more modest density dependence above jamming. A 
closer look to the numerical data indicates that the den¬ 
sity dependence of G is more pronounced above jamming 
than the one of B. Similarly to B , the sharp features 


of the shear modulus disappear rapidly when tempera¬ 
ture is increased above T = 10 -8 , and again the density 
dependence is very smooth when T > 10~ 6 . This in¬ 
dicates that the characteristic critical laws associated to 
the shear modulus near the jamming transition are easily 
smeared out by thermal fluctuations as well. 

The low-temperature crossover behaviour observed for 
G is again the signature of the zero-temperature critical¬ 
ity associated to the jamming transition. For ip < ipj and 
T —>• 0, the system explores the divergence of the shear 
modulus of Brownian hard spheres approaching jamming, 
which follows 


G ~ T(ipj - <p)~ K , ( 8 ) 

where k « 1.41 is a non-trivial critical exponent pH l47l . 
140 . On the other side of the jamming transition, jammed 
harmonic spheres lose shear rigidity as the jamming den¬ 
sity is approached from above G3, 

G~ {ip-ipj) 1/2 . (9) 

Although our numerical results are consistent with 
Eqs. @0, they are not precise enough to confirm that 
the exponent k is different from a previous estimate 
At = 3/2 Q, which is only marginally different from its 
recently predicted value k ~ 1.41. 


C. Ratio B/G of bulk to shear modulus: a 
signature of jamming criticality 

Whereas we noted that both B and G show qualita¬ 
tively similar sharp features in the vicinity of the jam¬ 
ming transition, we also stated that the precise values 
of the exponents characterizing their power law behav¬ 
ior are different. These differences stem from the fact 
that the behaviour of the bulk modulus B can be un¬ 
derstood from the evolution of the pressure, whereas the 
behavior of G is ruled by the evolution of the response 
of the system to a shear deformation. It is a specific sig¬ 
nature of the jamming transition that responses to shear 
and to compression differ maximally for isostatic pack¬ 
ings 0 , 1 , 0 , mm. 

Therefore, to clearly detect a quantitative sign of the 
jamming criticality, it is useful to analyse the behavior 
of the ratio B/G, which becomes infinite at the critical 
point. We combine our finite temperature data for B 
and G to follow the density dependence of B/G for vari¬ 
ous temperatures in Fig. 0 For our lowest temperature, 
T = 10 _s , we find that B/G is of order 5-10 far from 
the jamming density, but has a sharp maximum of or¬ 
der B/G « 100 when ip ss ipj. This behavior should be 
interpreted as a smooth version of the zero-temperature 
density dependence, which follows from Eqs. 0000 : 

B/G ~ {ipj - p) K ~ 2 , <p < ipj , (10) 

~ (¥>- <A7) _1/2 , V > <PJ, (11) 
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where k — 2 ss —0.59. Notice that the behaviour of B/G 
is now more symmetric around the jamming transition 
as the temperature prefactor disappears from the ratio 
B/G, but the critical exponents slightly differ on both 
sides of the transition (the divergence should be sharper 
for tp < ipj). Note also that the behavior of B/G is quan¬ 
titatively analogous to the behavior of the adimensional 
mean-squared displacement defined in Ref. jllj ]. There¬ 
fore, this figure demonstrates that the impact of thermal 
fluctuations on mechanical moduli and on single particle 
dynamics is actually identical, and mechanical moduli are 
in fact equally fragile against Brownian motion. 

We interpret the smoothened version of the symmetric 
divergence described by Eq. CD observed for T = 10 8 
in Fig. [T] as a ‘thermal vestige’ of the jamming transi¬ 
tion [18(, which is equivalent to the adimensional mean- 
squared displacement defined in Ref. [11]. These two 
quantities are both direct signatures of jamming criti¬ 
cality and should therefore be contrasted with a non¬ 
monotonic behaviour of the pair correlation functions 
mm which is instead a more general consequence of 
the particle softness and, as such, survives arbitrarily far 
from the critical point [5l|. Therefore we suggest that 
the observation of a large B/G ratio is a genuine sign 
that a particular material lies in the critical regime of 
the jamming transition, whereas a density-maximum in 
the pair correlation function is not. 

When the temperature is increased, the non-monotonic 
density dependence of B/G is rapidly erased by thermal 
fluctuations. For T = 10 -5 , we observe that B/G is 
nearly independent of the density and has a value of or¬ 
der 5 — 10 at all <p. These findings directly confirm that 
the jamming criticality is rapidly smeared out by thermal 
fluctuations. We also notice that even for very low tem¬ 
peratures, the range of densities where anomalous behav¬ 
ior associated to jamming can be observed is extremely 
narrow. These conclusions, obtained from the analysis of 
mechanical moduli, are in full agreement with previous 
conclusions drawn from the analysis of the mean-squared 
displacements Ell¬ 
in Ref. (2lj| . the bulk and shear moduli of a two- 
dimensional assembly of soft microgels were analysed, 
and their density evolution interpreted in terms of the 
power laws associated to the zero-temperature jamming 
criticality. Previous analysis of similar microgel systems 
has shown that these particles are quite soft, so that ther¬ 
mal fluctuations are of order T ss 10~ 6 —10 -4 , depending 
on the precise experimental system [THI, [52[ . The data re¬ 
ported in Ref. [2l| show that the ratio B/G is B/G ~ 3 
with a very weak density dependence. This is very much 
consistent with the physical interpretation that this sys¬ 
tem is far from being critical, which reinforces the gen¬ 
eral conclusion that very soft microgel systems are not 
good experimental systems to reveal thermal vestiges of 
the jamming transition. In Ref. [l5|, we suggested that 
emulsion droplets might be better suited for this task, as 
recently confirmed experimentally pT3? . 


IV. STATIC STRUCTURE FACTORS AND 
DIVERGING LENGTHSCALES 

In this section we define and study a number of static 
structure factors that can be probed in kinetically ar¬ 
rested colloidal materials in the presence of thermal fluc¬ 
tuations. From their low-wavevector analysis, we define 
lengthscales that diverge as the T = 0 jamming transi¬ 
tion is approached in the (T, ip) plane. 

A. Definitions of structure factors 

Because we know the position of each particle in each 
configuration, we can define a number of static structure 
factors from our particle packings at finite temperatures. 

The standard definition of the static structure factor 
is given by 

S ( k ) = ( 12 ) 

where we have defined the Fourier transform of the den¬ 
sity field as 

Pk = ^2exp(-ik ■ R 3 ), (13) 

j 

with Rj the position of particle j. We assume that all 
our packings are isotropic so that structure factors only 
depend upon wavevectors through their moduli k = |fc|. 

For a simple fluid at thermal equilibrium, the zero 
wavevector limit of S(k ) is directly related to the bulk 
modulus, and we have 

lim S(k) = ^- (Fluid). (14) 

This relation is only correct for liquid states [31]. We 
shall see below that a different analysis is needed for 
jammed solids. 

We consider systems that are lacking crystalline order 
but are nevertheless kinetically arrested. This implies 
that translational invariance is actually broken, and that 
we consider instead materials with long-range “amor¬ 
phous order” to use the language of glass theories (43lf54| . 
In practice, this means that the particle positions, and 
therefore, density fluctuations, can be naturally decom¬ 
posed into two different contributions, from which two 
distinct structure factors can be defined: 


S s(k ) = 

(15) 

s o(k) = 

(16) 

with 6p% = Pg — (Pk)- From the definitions (fl2l fl5l flGl) 
it is obvious that 

S(k) = S s (k) + S 0 (k), 

(17) 
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showing that we have obtained a decomposition of the 
structure factor in terms of a fluctuating part, Ss(k), 
and a configurational part, So(k). Physically, So(k) rep¬ 
resents the structure factor associated to the averaged 
position of the particles and is essentially independent of 
temperature and weakly dependent of density; it is an 
‘inherent’ property of the amorphous packing (55j . By 
contrast, Sg(k) represents the structure factor associated 
to the fluctuations of the particles away from their av¬ 
erage positions, and a strong temperature dependence is 
expected for this contribution, which should for instance 
vanish as T — > 0 for ip > ipj, when particles stop moving 
completely. 

When translational invariance is broken, as in crystals 
and glasses, the bulk modulus of the system is no longer 
related to S(k) as in Eq. (IhH) . but to the fluctuation part 
Sg(k) [3Cj. Within a conventional conventional elasticity 
theory where the elastic moduli are assumed to be inde¬ 
pendent of wavevector, an elastic body is characterized 
by longitudinal plane waves with the dispersion relation 

u> = k^J(B + | G)/p. These plane waves are thermally 
excited and follow the equipartition law, so that the fluc¬ 
tuating part of the density fluctuations is given by |(jj| 

Sg(k) = b (Continuum solid). (18) 

This is the fluctuation formula appropriate for connect¬ 
ing a static structure factor to the bulk modulus in a solid 
state. It is obviously distinct from the formula valid for 
fluids, Eq. m- Notice that the distinction between the 
two formula stems from the fact that translational in¬ 
variance is broken (in solids) or not (in fluids), but both 
formula rely on the fact that the system (fluid or solid) 
obeys the rules of equilibrium thermodynamics. 

It is also useful to provide a dynamic interpretation of 
the decomposition in Eq. m- Because the dynamics is 
arrested and the system only probes thermal fluctuations 
near a given metastable state, the dynamic structure fac¬ 
tor does not decay to zero at long times. The interme¬ 
diate scattering function is F(k,t ) = (p^(0)p_^(t))/N, 
so that F(k,t = 0) = S(k). In the long-time limit, 
density fluctuations are uncorrelated, and we directly 
find that F(k,t —> oo) = So(k), which is nothing but 
the collective Debye-Waller factor. Therefore, the fluc¬ 
tuation part of the static structure can be written as 
Sg(k) = S(k) — F(k,t —> oo), which quantifies the relax¬ 
ing part of the density fluctuations. 

In addition, we define two more structure factors as¬ 
sociated to the particle positions, which rely on the vec¬ 
torial character of the displacement field. In solid states, 
each particle vibrates around its average position. We 
can then define the displacement of each particle as 
Ui = Ri — (Ri), and the associated displacement field, 
expressed in the Fourier domain as 

U k = ^2ujexp(-ik- {Rj}). (19) 

3 


In the Fourier space, we can then decompose the displace¬ 
ment field into its longitudinal and transverse parts: 

% = ^ U L,k + %\fc’ (20) 

where u L £ = k ■ Ur, and k = k/\k\ is the unit vector in 

the direction of k. Using these fields, we can finally define 
the longitudinal and transverse correlation functions 

S L (k) = 

Sr(k) = — {u t £ ■ u T _j~). (21) 

Following again the approach of conventional elasticity 
theory, we now have the following expressions j30| 

SM = ( 22 ) 

St (A;) = (Continuum solid). 

These expressions directly show that we can calculate 
the shear modulus G from the low-wavector limit of 
Sr(k) 0- It is the approach that has been employed to 
obtain the data shown in Fig. Q] in Sec. IIIII 

Before using these definitions we wish to comment that 
the various structure factors defined in this section cor¬ 
respond to different mathematical ways of decomposing 
the particle positions into an average and a fluctuating 
contributions. In the first approach, we performed the 
decomposition in the Fourier domain, whereas the second 
approach deals with position space. It should therefore 
come as no surprise that the longitudinal part Si,{k) can 
be related to the fluctuating part of the static structure 
factor Sg{k). In terms of the particle displacements Uj, 
the fluctuating part of the density field can be expressed 
as 

5pj: = ^^(—ik ■ Uj + 0{k 2 )) exp (—ik ■ (Rj))- (23) 

3 

Comparing this expression to Eq. (HU), we get Spj: ~ 
—iku L t at the lowest order in k. Therefore we conclude 
that at this order 

S L (k) « S s (k), (24) 

and both approaches actually carry equivalent physical 
content, as they should. By a similar reasoning, we find 
that 

(Pk) =J2^p(-ik-(Rj))(l-k 2 ((k-Uj) 2 )/2+---), (25) 

3 

which shows that an accurate estimate of the configura¬ 
tional part So(k) can be obtained by computing the struc¬ 
ture factor of the average particle positions, a method 
that could prove convenient for experiments using parti¬ 
cle imaging. 
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FIG. 2: (Color online) (a) Fluctuation part of the static struc¬ 
ture factor Ss(k) at T = 10 -8 for various densities, (b) Same 
data with vertical axis scaled with the factor (B + | G)/(pT) 
expected from the usual plane-wave behavior, (c) The hor¬ 
izontal axis is also rescaled to obtain the best data collapse 
and extract the longitudinal length scale £l, Eq. m- 


B. Longitudinal fluctuations 

We first discuss the behavior of Sg(k) in connection 
with the bulk modulus. In practice, we have calcu¬ 
lated Ss(k) in the following way. We first calculate 
the density field p g for each instantaneous configurations 
of particles, at the lattice points of the first Brillouin 
zone, k = ( k Xl k yi k z ) in k -space where k a = n a Tt/L for 


a = x,y,z and n a an integer. We then perform a time 
average to obtain (p~) for each Fourier component. Us¬ 
ing pj: and ( p j) we calculate Spj: and obtain Sg(k) in a 
straightforward manner. We finally perform a circular 
average to obtain the desired Sg{k). 

We plot the numerical results for Sg(k) at T = 10 ~ 8 
and various densities in Fig. [2] The overall amplitude of 
Ss(k) strongly decreases when the system is compressed. 
This is expected, because particles move less and less 
when density is increased [ll], and the overall amplitude 
of the fluctuations gets smaller. At all densities, we also 
find that Sg(k) has a well-defined k —> 0 limit, and that it 
increases strongly with k , with a well-defined first diffrac¬ 
tion peak corresponding to the interparticle distance at 
k ~ k o = 27r/cr, which reflects the liquid-like structure of 
amorphous packings at the particle scale. 

Two useful informations are contained in these struc¬ 
ture factors. We first concentrate on the k —0 limits, 
where the relation Eq. (1181) derived from continuum the¬ 
ory is expected to become valid, 

iTtA^'sTfc' (26) 

To verify this fluctuation formula, we use this expression 
and replot the same data in a scaled form, observing the 
fc-dependence of the quantity Sg(k)(B + | G)/(pT), as 
shown in Fig. [2] (middle). The bulk and shear moduli 
are obtained from independent measurements, shown in 
Fig. CD Note that since our systems are characterized by 
large B/G ratio, as discussed in Sec. IIIII it means that 
the bulk modulus gives the major contribution to the 
(B + |G) factors in all these expresssions, and the term 
G is almost always negligible. 

The numerical results show that the quantity 
Sg(k)(B + | G)/(pT) clearly approaches unity as fc —> 0, 
as it should. This observation implies that standard 
equilibrium relations between the mechanical moduli and 
static structure factors are valid at large lengthscales for 
jammed packings, provided the appropriate fluctuating 
part of the density fields are analysed, instead of the to¬ 
tal S(k). 

The rescaled plot shows moreover that Sg(k) not only 
contains useful information at k —> 0, but that its finite 
k behaviour is also relevant. The conventional elastic¬ 
ity expression in Eq. m does not provide any useful 
fc-dependence. In other words, when all the longitudinal 
waves are the plane waves predicted by conventional elas¬ 
ticity theory, the quantity Sg(k)(B + | G)/{pT) should 
be unity. Thus, deviations from unity can be interpreted 
as an indicator for the breakdown of the description by 
the usual plane wave with the above dispersion relation. 
Interestingly we find that the rescaled data in Fig. [2] 
(middle) show fmite-fc deviations that depend strongly 
on the density, and are maximal at the jamming density, 
so that the deviations from conventional elasticity have 
a remarkable non-monotonic density dependence. Inter¬ 
estingly, near the jamming density, ip « ipj, we observe 
a very clear power law behaviour, Sg(k) oc k 2 . 
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To characterize quantitatively these deviations and 
their apparent relation with the jamming criticality, we 
propose the following scaling analysis of these data. The 
above description suggests the existence of a non-trivial 
correlation lengthscale, separating two distinct be¬ 
haviours: Sg{k^L 1 ){B + | G)/(pT) ~ 1, and Sg(l 
k^L fco£i)(-B + | G)/(pT) oc k 2 . Therefore, we de¬ 
termined numerically the longitudinal correlation length- 
scale assuming the scaling form 

Ss(k) « (27) 

B+ 3 G 

where F(x) is a scaling function of the form F(x C 1) = 
1 and F(x 1) oc x 2 . Physically, this expression implies 
that is a characteristic length scale below which the 
conventional elasticity description of longitudinal parti¬ 
cle displacements breaks down. A diverging correlation 
length implies that the usual plane wave description 
does not apply at any length scale. 

The results of this scaling analysis are shown in Fig. [2] 
(bottom). The data collapse is acceptable, but it is diffi¬ 
cult to confirm its validity, because the obtained length- 
scale is quite large, and a larger system size would be 
required for a more accurate determination of this quan¬ 
tity, especially close to <pj at very low temperatures. The 
evolution of the obtained longitudinal lengthscale with 
ip and T is discussed below in Sec. IIVDI 


C. Transverse fluctuations 

We now repeat the analysis of Sec. lIVBl for the evolu¬ 
tion of Srik) in connection with the shear modulus. We 
plot Sr(k) at T = 10 -8 and various densities in Fig. [3] 
(top). As for Ss(k), we find that the overall amplitude of 
Sr{k) decreases rapidly with the density, with an overall 
wavevector dependence similar to the one of Ss(k). 

We now use the fluctuation formula for the transverse 
fluctuations (i3 | 

lim Sr(fc) = ^, (28) 

k-¥ 0 Lx 

which is the method we have employed to measure the 
shear modulus presented in Fig. [T| Therefore, by con¬ 
struction, when we replot the quantity St(A:)(G/ pT) 
in Fig. [2] (middle), the data extrapolate to unity when 
k —> 0. 

This figure shows again that the deviations from con¬ 
ventional elasticity have a striking non-monotonic den¬ 
sity dependence, and that deviations are maximal when 
ip is close to tpj. Following the analysis of longitudinal 
fluctuations, we again hypothesize a scaling behavior for 

Sr(k): 

S T (k) « ^ff(fcfr), (29) 

where H{x) is a scaling function of the form H(x <C 
1) = 1 and H{x 1) oc x 2 . We show in Fig. [3] (bottom) 



FIG. 3: (Color online) (a) Transverse part of the displacement 
structure factor Sr(k) at T = 10~ 8 for various densities, (b) 
Same data with vertical axis scaled with the factor G/(pT) 
expected from the usual plane-wave behavior, (c) The hor¬ 
izontal axis is also rescaled to obtain the best data collapse 
and extract the transverse length scale £t, Eq. ini). 


a collapse of the numerical data according to Eq. HMD 
which allows us to determine numerically a lengthscale 
£t which delimits the validity of the usual plane wave 
description of transverse fluctuations of the particle dis¬ 
placements. This critical scaling law implies again that 
the usual plane wave description does not apply at any 
length scale when £t diverges, which is expected exactly 
at the jamming transition. 
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FIG. 4: (Color online) Density dependence of £l (a) and 
£t (b) for various temperatures. The near-critical non¬ 
monotonic density dependence observed for T = 10~ 8 be¬ 
comes a smooth variation for T = 10 -5 when jamming crit¬ 
icality is smeared by thermal fluctuations. The usual plane 
wave description of vibrational motion is therefore excellent 
away from criticality, T > 10~ 6 and/or |ip — pj | > 0.02. 


D. Transverse and longitudinal lengthscales 

We have performed the scaling analysis outlined in 
Secs. IIV Bl and ITVCl for different temperatures and have 
been able to extract the density and temperature depen¬ 
dences of the characteristic length scales £l and £r. The 
results are presented in Fig. [I] 

Both lengthscales behave qualitatively similarly. At 
T = 10 -8 , they strongly depend on the density: they 
have a clear maximum for p s=s pj and become of order 
unity far away from the jamming transition. In particu¬ 
lar, we observe that £l becomes comparable to the sys¬ 
tem size (L = 36.2) for T = 10~ 8 , which explains why 
its numerical determination is difficult. By contrast, the 
maximum value reached by £t is more modest (or order 
£t ~ 10), explaining why the data collapse for Sr(k) is 
more convincing than the one for S$(k). 

When the temperature is increased, the density max¬ 


imum observed for the correlation lengthscales is much 
less pronounced, and nearly disappears when T = ICR 5 
where £^ and £t have uninteresting density dependences 
and remain microscopic, £l,t ~ 2 — 5. The conclusion 
is that in this regime, which is relevant for a number of 
experimental situations for colloidal systems, continuum 
mechanics actually represents an excellent approxima¬ 
tion down to microscopic length scales. In other words, 
“anomalous” or “soft” modes, which exist over arbitrary 
length and frequency scales at the jamming transition 
where correlation lengthscales and timescales are infinite, 
are strongly suppressed by moving away from jamming 
criticality. 

We now compare our measurements of the lengthscales 
£l and £t to similar lengthscales measured earlier in the 
literature. A first relevant comparison is with the re¬ 
sults in Refs. m, where a characteristic length scale 
for longitudinal and transverse plane waves at a specific 
frequency w* = u>*(p) for p > pj and T = 0 were mea¬ 
sured. In this approach, w* is a characteristic frequency 
where the vibrational density of state exhibits anomalous 
(nearly frequency-independent) behavior 13, H- The ob¬ 
tained longitudinal and transverse length scales £2 and 
£2 measured from this protocol are predicted to diverge 
as £2 ( ( P~ 1 Pj)~ 0 ' 5 and £2 °c (v? — <Pj) - 0 ' 25 (Ell, and the 

latter behavior is directly confirmed by numerical mea¬ 
surements (data for £2 were not shown). 

Although it is tempting to assume that £l and £t have 
similar physical content as £2 and £2, power law fits to 
our measurements yield exponents that are not consis¬ 
tent with the predicted 0.5 and 0.25 (we find that 0.7 
and 0.5 fit our data better). However, the lengthscales 
observed in our measurement is so large that their pre¬ 
cise determination is challenging. Much larger systems 
are needed to fully settle this issue. Furthermore, our 
lack of knowledge about the precise form of the scaling 
functions F(x) and H(x) may also affect the quality of 
our estimates for these lengthscales. We wish to raise 
the possibity that the two sets of lengthscales are not 
fully equivalent because we directly defined characteris¬ 
tic lengthscales where the usual plane wave descriptions 
stop working, whereas Silbert et al. focused on a specific, 
density-dependent frequency u>* [Mill- 

In a more recent work, Wang et al. [Hl| have deter¬ 
mined numerically the so-called Ioffe-Regel frequencies 
and lengthscales. These are respectively defined as the 
timescales and lengthscales characterizing the disappear¬ 
ance of plane waves in the collective dynamic structure 
factors. We may expect that the lengthscales £l and £t 
that we have defined above carry a similar physical con¬ 
tent to the Ioffe-Regel lengthscales II and It analysed 
in Ref. [58;] • Although some of the numerical results of 
Wang et al. are consistent with ours, their final conclu¬ 
sions differ qualitatively from ours. In particular, their 
analysis suggests that in the low-temperature limit It di¬ 
verges at <p = pj and transverse plane waves do not exist 
in the hard sphere regime p < pj and T —» 0. Instead, 
we observe that £t becomes finite on both sides of the 
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FIG. 5: (Color online) Configurational part of the static struc¬ 
ture factor So(k) at p = 0.652 > pj and several temperatures. 
Circles, squares, and upper and lower triangles are represent 
the data obtained at finite temperatures for N = 64, 000, di¬ 
amonds indicate data obtained at T = 0 for N = 512, 000. 
Solid line indicates the hyperuniform behaviour of density 
fluctuations at large scale, So(k) oc k. The hyperuniformity 
appears essentially independent of temperature. The inset is 
the zoom out of the main panel. 


jamming transition. We suspect that their interpreta¬ 
tion is incorrect and stems from extrapolating numerical 
measurements performed at <p > ipj to the hard sphere 
glass. Our results show instead that longitudinal and 
transverse vibrations do exist in the hard sphere regime, 
and the associated lengthscales and timescales actually 
become microscopic as the density is decreased away from 
<pj. Contrary to the claim in Ref. 58}, the hard sphere 
glass is not qualitatively different from other amorphous 
materials. 


V. HYPERUNIFORMITY OF THE 
CONFIGURATIONAL STRUCTURE FACTOR 


In this section, we analyse the structure factor asso¬ 
ciated to averaged particle positions and show that this 
configurational part So(k) reveals a nearly-hyperuniform 
behaviour at large length scales. We also show how the 
analysis of So(k) allows us to elucidate the physical con¬ 
tent of the related nonequilibrium index and nonequilib¬ 
rium lengthscale measured in amorphous materials. 


A. Hyperuniformity is unrelated to jamming 
criticality 

In Sec. II V Al we decomposed the structure factor S(k) 
into a fluctuating part, analysed in Sec. IIVBI and a 
configurational part S'o(fc), which is the subject of the 
present section. 


The first question we ask is whether the configurational 
part So(k) is sensitive to the temperature. In Fig. [5j we 
show numerical results for So(k) at a constant density 
p = 0.652 (slightly above pj) at several temperatures 
from T = 10“ 8 up to T = 10 -5 . The inset shows that 
S'o(fc) has the usual wavevector dependence of a liquid 
structure factor, with a broad first diffraction peak near 
k « k 0 . Clearly, the temperature dependence appears 
negligible in this representation. 

In the main panel, we zoom on the low-wavevector 
behaviour in order to reveal a possible effect of the ther¬ 
mal fluctuations. Within the accuracy of these computa¬ 
tions, we find again no visible temperature dependence 
for So(fc). This clearly confirms that thermal fluctuations 
mainly contribute to the fluctuating part of the density 
fluctuations, whereas the averaged component of density 
fluctuations is essentially unaffected by the temperature, 
at least for the range of wavevectors shown in Fig. [5] 
We expect more changes to occur at very large wavevec¬ 
tors, k k 0 , where the sharp features associated to the 
pair correlation function at contact produce long-ranged 
oscillations (59} • 

An interesting behaviour is observed for the low-fc be¬ 
haviour of So(k) in Fig. [5] In the linear scale chosen 
for representating these numerical measurements, we ob¬ 
serve that So(k) oc k, for k < 2. This “anomalous” lin¬ 
ear behavior with an apparent vanishing of So(k —> 0) 
has been termed hyperuniformity |22l [23} . Hyperuniform 
density fluctuations have been reported in simulations 
of sphere pac kings a t the jamming transition both nu¬ 
merically [22L l24i-[26l and in experimentally-constructed 
granular packings j25|. In colloidal systems, signs of hy¬ 
peruniformity are much weaker [27H29j |. 

The important conclusion that we can draw from the 
absence of temperature dependence in the data shown 
in Fig. [5] is that hyperuniformity appears extremely ro¬ 
bust against thermal fluctuations, and can in fact eas¬ 
ily be observed even for our highest studied tempera¬ 
ture T = 10~ 5 . This observation is in striking contrast 
with all other observations reported earlier in this pa¬ 
per related to the jamming criticality. Whereas quanti¬ 
ties such as mechanical moduli and correlation length- 
scales are rapidly smeared out by thermal fluctuations, 
hyperuniformity appears rather insensitive to tempera¬ 
ture. This strongly suggests that jamming criticality and 
hyperuniformity are unrelated concepts and have distinct 
physical origins. 

We emphasize that the decomposition of the structure 
factor S(k) as the sum of two terms So(k) and Ss(k) 
indicates that the total structure factor is related to 
the isothermal compressibility only through Sg(k) whose 
wavevector dependence shows no anomalous dependence, 
see Fig. [2] On the other hand, we find that So(k) is 
characterized by a hyperuniform linear behaviour at low 
wavevector, but this configurational contribution is un¬ 
related to the compressibility. Thus, we conclude that 
hyperuniformity (related to So) cannot be a logical con¬ 
sequence of a vanishing compressibility (related to Sg ) of 
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the packing. Of course, to observe hyperuniformity in 
the total structure factor S(k) rather than in the config¬ 
urational part So(k), it is necessary that Sg(k ) -C So(k), 
which happens when the compressibility becomes very 
small. Working close to the jamming transition is there¬ 
fore a practical rather than a fundamental issue , which 
is no longer needed when working directly with So(k). 
In that sense, the discovery that some amorphous ma¬ 
terials become liyperunifornr very close to the jamming 
transition may appear coincidental [22j. 

These findings are experimentally relevant, because 
they mean that hyperuniformity, unlike jamming critical¬ 
ity, can actually be observed across a large region of the 
(T, p) phase diagram for soft colloids. Quite surprisingly, 
the experimental literature seems to suggest exactly the 
opposite since a number of experiments have reported 
signatures of jamming criticality in soft colloids [H, HU 
(which, we argue, are actually taken far from critical¬ 
ity), whereas only weak signs of hyperuniformity have 
been reported 27H29| (which, we argue, should be easily 
observable). This might be due to the difficult experi¬ 
mental constraint that measuring S'o(fc) at low k requires 
detecting the position of a large number of particles with 
a large precision. 



FIG. 6: (Color online) Evolution of So(k) with the density 
above the jamming transition in a log-log representation. 
These data are obtained by averaging over 8 independent 
packings with N = 512, 000 particles. A clear hyperuniform 
linear ^-dependence (shown as a full line) is obtained over 
a broad range of wavevectors when tp < 0.7 with only weak 
density dependence, but the data saturate to a finite value as 
k ->• 0. 


B. Density dependence and deviations from strict 
hyperuniformity 

Having established that temperature does not affect 
much the observation of hyperuniformity, we then dis¬ 
cuss the effect of the density by analysing data obtained 
directly at T = 0. This approach is useful, because it al¬ 
lows us to disentangle temperature and density effects. In 
addition, we can study at T = 0 much larger systems in 
order to analyse whether hyperuniform behaviour can be 
observed over arbitrarily large lengthscales. To this end, 
we employ a distinct measurement technique of S’o(fc) at 
T = 0 with a larger number of particles, N = 512, 000, 
and use 8 independent packings to reduce the statisti¬ 
cal error. Our numerical methodology was described in 

Sec. ITT] 

Using this larger systems at T = 0, we confirm in 
Fig. [5] that So(k) for this density is in excellent agree¬ 
ment with the finite temperature results obtained with 
the smaller packings, although of course the statistical er¬ 
ror is greatly reduced. The agreement between these two 
independent sets of data confirms that S’o(fc) is largely 
independent of temperature. 

Having shown that temperature plays no role, we can 
now analyse the density dependence of these observa¬ 
tions. In Fig. [6] we plot So(k) measured at T = 0 at 
various densities between tp = 0.80 much above jam¬ 
ming, down to tp = 0.652 just above pj. We now use 
a log-log representation of the results. At p = 0.80, 
So(k) is almost constant over a large range of wavevec¬ 
tors, So(k) « 4 • 10 -3 for k < 1 , and a hyperuni¬ 
form behaviour cannot be observed. This behaviour of 


over-compressed packings is consistent with observations 
made in other glass-forming materials, such as simple 
Lennard-Jones glasses, where hyperuniformity is not ob¬ 
served either [261] . 

For p = 0.7, a linear behavior, S'o(fc) « 0.004A:, already 
appears in wide k region, even though \p — pj\ ss 0.055 
is still quite large. More surprising is the observation 
that the data for p = 0.652, 0.66 and 0.68 (respec¬ 
tively corresponding to |p — pj\ « 0.0063, 0.015, and 
0.034) are essentially the same, and are characterized 
by a broad range of wavevectors with linear dependence, 
So(k) ~ 0.004fc, although the data saturate at very low 
A: to a finite value So{k —> 0) « 1.4 • 10 -3 . These re¬ 
sults indicate that hyperuniformity is a robust feature of 
S 0 (k), in the sense that it is weakly dependent on the 
density and does not require fine-tuning the volume frac¬ 
tion to the jamming density pj, confirming that the two 
concepts are distinct. 

However, it should also be noted that a strict hyperuni¬ 
formity S(k) oc k can not be observed down to arbitrarily 
small wavevectors, and deviations appear below k ~ 0.4, 
which corresponds to a large lengthscale ~ 15cr. This 
surprising saturation effect has not been reported before, 
although we notice that previous literature [22|, H(| in¬ 
dicates that the smallest So(k —»• 0) values achieved in 
computer simulations are always of the order 10 -3 or 
more, which is consistent with our own results. This sat¬ 
uration would not be observed if the data in Fig. [0] where 
plotted in a linear scale. 

Our analysis shows that this saturation effect is clearly 
not due to thermal fluctuations (we work at T = 0). 
This does not stem from sample-to-sample fluctuations 
either, because all 8 samples show a saturation of simi- 
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lar amplitude. Finally, the saturation does not seem to 
depend on density, at least for p > ipj. We cannot ac¬ 
cess the regime <p < ipj using energy minimization, but 
we note that a marked density dependence was reported 
for the structure factor of hard spheres approaching the 
jamming density from below in Ref. f22]. However, the 
total structure factor S(k) was measured in that study, 
which contains a density-dependent contribution associ¬ 
ated to the fluctuating part Ss(k ), which vanishes as ipj 
is approached. It would therefore be very interesting to 
measure directly So(k) in the hard sphere glass for very 
large system sizes. We have performed exploratory sim¬ 
ulations with moderate system sizes in this regime and 
find a weak density dependence of So(k ) when p > 0.62, 
but larger systems are needed to analyse more finely the 
behaviour at very small k. 

Finally, we remark that it is difficult to provide a phys¬ 
ical explanation for the existence of the observed devia¬ 
tions from strict hyperuniformity, mainly because there 
is no deep physical reason to expect perfect hyperunifor¬ 
mity in these systems in the first place. Among possible 
factors that could be investigated are the role of a fi¬ 
nite density of rattlers in the packings and the role of 
the specific protocol that is used to prepare the particle 
packings, which both could influence the measured value 
of Sq( k —> 0). Such studies are beyond the scope of the 
present work. 


C. Analysis of the “nonequilibrium index” 


We showed in Sec. II V Bi t. hat the isothermal compress¬ 
ibility is directly related, for solids at thermal equilib¬ 
rium, to the fluctuating part of the structure factor via 
Eq. (fl8l) . A direct consequence is that the compressibil¬ 
ity is thus not related to the total structure factor S(k ) 
via the relation valid for equilibrium fluids , Eq. (1141) . To 
quantify the difference between fluid and solid states, the 
concept of a “nonequilibrium index”, A, was introduced 
and studied both for hard sphere glasses [321 ] and for other 
types of amorphous materials [33J ■ 

We now show that the decomposition provided above 
for the structure factor allows us to elucidate the phys¬ 
ical content of X. Using the notations introduced in 
the present work, the nonequilibrium index X is defined 
as 321 


x = , im ?m _ ! 

fc->o pT 


(30) 


By construction, X = 0 for a fluid at thermal equib- 
rium, see Eq. m- Since A' is defined as the ratio 
between fluctuations and response functions, its func¬ 
tional form is also reminiscent of the effective temper¬ 
ature and fluctuation-dissipation ratio that characterize 
the nonequilibrium of aging and driven glasses [34j . The 
main difference between the two types of quantities is 
that A' refers to static fluctuations, whereas effective 
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FIG. 7: (Color online) Temperature dependence of the 
nonequilibrium index A' measured in two model glasses (sym¬ 
bols). Full lines are from Eq. (1311) . the prediction obtained by 
assuming that equilibrium fluctuation relations are satisfied 
for the glass. 


temperatures are defined from time-dependent correla¬ 
tion and response functions. The non-zero value mea¬ 
sured for X in hard sphere glasses was interpreted as 
a demonstration that the 11 jammed glassy state is funda¬ 
mentally nonequilibrium in nature" [32[. The simulations 
indicate that X grows rapidly when hard sphere glasses 
are compressed towards ipj, or when the temperature is 
decreased below the glass transition temperature T g of 
model glass-forming systems, such as Lennard-Jones and 
Dzugutov glasses 133 ] . 

The decomposition of the structure factor proposed in 
Eq. (1171) provides us with two important informations. 
First, we have shown that equilibrium fluctuation rela¬ 
tions are perfectly obeyed in the solid phase for static 
quantities. This result implies that the nonequilibrium 
nature of the glass cannot be revealed by a fluctuation 
formula based on static density fluctuations and sug¬ 
gests, in fact, that the introduction of a “nonequilib¬ 
rium” index to characterize static density fluctuations 
is unnecessary. This conclusion is in qualitative agree¬ 
ment with the two-temperature scenario for the nonequi¬ 
librium dynamics of glasses, where short-time fluctua¬ 
tions and response are typically found to obey equilib¬ 
rium fluctuation-dissipation relations [34, 35]. 

Second, the combination of Eqs. mm provides pre¬ 
dictions for the leading behaviour of the nonequilibrium 
index in various systems. For glass-forming models with 
continuous interactions, we can assume that So(k —> 0) 
and the bulk modulus are weakly temperature depen¬ 
dent deep in the glass phase |59|, so that in the low- 
temperature limit, one gets 


X(T « T g ) 


So(k -A 0)B 1 

flf “ T' 


(31) 


In Fig. [7] we confirm that the low-temperature behaviour 
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FIG. 8: (Color online) Temperature dependence of the 
nonequilibrium lengthscale £ neq measured in two model 
glasses (symbols). Full lines are from Eq. (13-41) . the prediction 
obtained by assuming the equilibrium fluctuation relations 
are satisfied for the glass. 


of the nonequilibrium index measured numerically in 
Ref. 0 is consistent with our prediction in Eq. m 
that it should diverge as 1/T as T —>■ 0. This leading 
temperature behaviour stems from the fact that S(k) in 
the definition of X in Eq. (HU contains a ‘frozen’ contri¬ 
bution, So, which does not vanish at T = 0. 

For hard sphere glasses, the leading asymptotic be¬ 
haviour of the nonequilibrium index depends strongly of 
the hypothesis made regarding the behaviour of So(fc) 
very close to ipj. Assuming that hyperuniformity is only 
approximate, as we report in Fig. [G] one would then pre¬ 
dict that X in Eq. (l30l) is dominated by the divergence of 
the bulk modulus, yielding X m (ipj — ip)~ 2 . In Ref. (32] 
a linear decrease, So(k -4 0) ~ (ip — <pj), was assumed 
for So, which turns into a different divergent behaviour, 
X ~ (ipj — </?) , for the nonequilibrium index. Nu¬ 

merically, we expect that X exhibits a crossover between 
these two power law regimes as ipj is approached, which 
could be difficult to analyse. 


D. Analysis of the “nonequilibrium lengthscale” 


We finally discuss the concept of a nonequilibrium 
lengthscale, £ neq , defined again for amorphous materi¬ 
als from the behaviour of the static structure structure. 
The nonequilibrium lengthscale is defined as 


£neq = [—c(fe -A 0)] 


1/d 


[pS(k -a 0)] 


-1/d 


(32) 


where c(k) is the direct correlation function IdTl and d 
is the dimensionality of the system. In the final part of 
Eq. (1521) . we have assumed that S(k — > 0) <C 1. In com¬ 
puter simulations, it is found that £ neq grows as tempera¬ 
ture is decreased below T g in model glasses [33:] and satu¬ 
rates to a finite value as T -A 0, whereas it is predicted to 
diverge as ip —»• ipj in hard sphere glasses 0. As such, 
it is interpreted as a growing static length scale that is 
potentially relevant to characterize the structure of the 
glassy state. In this view, hard sphere glasses would 
therefore be somewhat “special” since they would have 
a diverging static lengthscale, whereas glasses with con¬ 
tinuous interactions would exhibit a non-diverging static 
lengthscale. 

Because this lengthscale directly follows from the low- 
k behaviour of the structure factor, our decomposition 
(da into two distinct contributions is again relevant to 
understand the physical content of the nonequilibrium 
lengthscale, which we can rewrite 


£neq — 


1 -1/d 


pS 0 (k -a 0) 


B + lG 


T 


(33) 


For Lennard-Jones and Dzugutov models, we again ex¬ 
pect that So(k -A 0) and the mechanical moduli are 
weakly dependent on temperature far below T ff , so that 
the leading temperature dependence of the nonequilib- 
rium lengthscale is transparent in Eq. (l33l) . and should 
be of the form 


£neq « (a + bT)- 1 ^, (34) 

where a and b are some constants. In Fig. [8] we confirm 
that this prediction describes the numerical data very 
well for two glass-formers, showing that the growth of the 
nonequilibrium lengthscale with decreasing temperature 
can in fact be fully understood by assuming that den¬ 
sity fluctuations obey equilibrium fluctuation formula. 
In essence, therefore, the growth of the nonequilibrium 
lengthscale in the glass phase reflects the competition 
between the configurational and fluctuating parts of the 
static structure factor, which have different temperature 
dependences: the former is essentially constant and re¬ 
flects the ‘inherent’ structure of the glass, the latter stem¬ 
ming from vibrational motion and is thus proportional to 
temperature, as captured in Eq. (1M1) . 

For hard sphere glasses, the behaviour of the nonequi¬ 
librium length would again depend sensitively on the be¬ 
haviour of So(k) near ipj. Assuming that hyperunifor¬ 
mity is only approximate, the nonequilibrium lengthscale 
would grow strongly as the glass phase is entered and the 
compressibility decreases, but its growth would saturate 
to a value £ neq « [pSo(k -A 0)] _1 / d « 8.2, using nu¬ 
merical values from Fig. 0 Interestingly, this saturation 
value is close to the value £ neq (T —> 0) ~ 7.5 found for 
the three-dimensional Dzugutov glass-former in Fig. [8] 
which could support the idea that hard sphere glasses are 
not a “special” type of glass-former. For a strictly hy¬ 
peruniform hard sphere system, on the other hand, the 
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nonequilibrium lengthscale would diverge as ip —> ipj, as 
predicted in Ref. j32l ]. 

VI. CONCLUSION 

In this work, we have analyzed the density and tem¬ 
perature dependences of mechanical moduli and several 
types of structure factors in a model system of soft har¬ 
monic spheres in the vicinity of the jamming transition. 

We have shown that thermal fluctuations very quickly 
erase several signatures of the criticality associated to 
jamming, in agreement with earlier work related to sin¬ 
gle particle dynamics & We showed that the bulk 
modulus, the shear modulus, the longitudinal and trans¬ 
verse lengthscales rapidly acquire a ‘normal’ behaviour 
typical of ordinary solids, whereas the large lengthscales 
and timescales associated to the isostatic jamming crit¬ 
ical point are only observed in a narrow region of the 
(T, <p) phase diagram. We conclude that most colloidal 
experiments to date have hardly been able to probe the 
jamming criticality, nor have the thermal vestiges of the 
jamming transition that result from the existence of non- 
microscopic lengthscales and timescales been observed. 
These conclusions suggest that the soft and hard sphere 
glasses that are commonly studied experimentally essen¬ 
tially behave as ordinary solids where the usual plane 
wave description holds down to small length scales, as 
concluded from a very recent experimental study 16 Of . 
Therefore, we hope that our results will encourage further 
experimental investigations of these issues, for instance 
using emulsion droplets [53j , or core-shell microgel parti¬ 
cles |61] . 

A second major finding in our study is that density 
fluctuations for jammed colloidal systems follow the laws 
of equilibrium thermodynamics and their study does not 
reveal the nonequilibrium nature of glasses. Our anal¬ 
ysis is based on a decomposition of density fluctuations 
in a configurational and fluctuating parts. Whereas the 
fluctuating part is directly related to mechanical moduli 
via equilibrium fluctuation formula, we found that the 
configurational part is essentially independent of both 


density and temperature in a rather broad range of pa¬ 
rameters. The decomposition into these two compo¬ 
nents allows us to elucidate the behaviour reported in 
earlier numerical studies for the nonequilibrium index 
and nonequilibrium lengthscales characterizing amor¬ 
phous materials. Based on these observations, we have 
suggested that hyperuniformity observed in the configu¬ 
rational structure factor is unrelated to the compressibil¬ 
ity, and therefore to the jamming criticality. 

These results raise some interesting questions. It has 
been established numerically that the same jamming crit¬ 
icality is observed for packings with very different prepa¬ 
ration protocols [12]. Our observation that a strict hy¬ 
peruniformity is not observed in our packings suggests 
that the value of So(k —» 0) could very well be affected 
by the nonequilibrium protocol used to prepare packings. 
One could for instance hypothesize that a packing pre¬ 
pared with a slower annealing could be more hyperuni¬ 
form than one produced via brutal compression. This 
raises the appealing possibility that the nonequilibrium 
lengthscale £„eq measured either at T = 0 (for continuous 
potentials) or at infinite pressure (for hard spheres) truly 
encodes some non-trivial information about the glassy 
state [32] • If correct, it would mean that it is not really 
the temperature or density dependences of £ neq which 
truly matter, but rather its evolution for different prepa¬ 
ration histories. Therefore, we believe that it would be 
interesting to understand better the physical content of 
this quantity in various glassy materials prepared using 
various thermal histories. 
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